#######################################
#########FACTORIAL EXPERIMENT##########
#######################################


#LOAD /INSTALL PACKAGES
pacman::p_load(
  ggplot2, gridExtra, cregg, margins, prediction, foreign, binom, Barnard, haven, Rmisc, sjmisc, magrittr, tidyverse, jtools, MASS,
  ggsci, 
  patchwork
)

####ANALYSIS WITH SPEEDERS####

### SET WORKING DIRECTORY TO FILE LOCATION
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#LOAD DATA
paper_risp_f<-read.dta("immigration_speeders.dta")

#EXPERIMENTAL FACTORS
paper_risp_f$approved_asy  <- factor(paper_risp_f$approved_asy)
paper_risp_f$exp_1  <- factor(paper_risp_f$exp_1)
paper_risp_f$exp_2  <- factor(paper_risp_f$exp_2)
paper_risp_f$exp_3  <- factor(paper_risp_f$exp_3)
paper_risp_f$ideo  <- factor(paper_risp_f$ideo)
paper_risp_f$ideo2  <- factor(paper_risp_f$ideo2)

#GLM MAIN EFFECTS (NOT WEIGHTED)
m3way <- glm(approved_asy ~ exp_1*exp_2*exp_3, family=binomial(), data = paper_risp_f)
summary(m3way)

###PREDICTIONS MAIN EFFECTS
prediction(m3way, at = list(exp_1 = c("Low skilled", "Skilled"))) %>% 
  summary -> main_skilled

prediction(m3way, at = list(exp_2 = c("Looking for job", "Fleeing from war"))) %>% 
  summary -> main_reason

prediction(m3way, at = list(exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> main_origin


### GGPLOT MAIN EFFECTS 
ggplot(data = main_skilled,
       aes(x = `at(exp_1)`,
           y = Prediction,
           shape = `at(exp_1)`,
           color = `at(exp_1)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
                   axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot1

ggplot(data = main_reason,
       aes(x = `at(exp_2)`,
           y = Prediction,
           shape = `at(exp_2)`,
           color = `at(exp_2)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("seagreen4", "palevioletred4")) +
  scale_shape_manual(name="", values = c(15, 18)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
                   axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot2

ggplot(data = main_origin,
       aes(x = `at(exp_3)`,
           y = Prediction,
           shape = `at(exp_3)`,
           color = `at(exp_3)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("steelblue4", "antiquewhite4")) +
  scale_shape_manual(name="", values = c(4, 8)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
                   axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot3

plot1 + plot2 + plot3 

ggsave("figureApp4.png", dpi = 600)


#####MODELS SPILL-OVER SPEEDERS#####

### SET WORKING DIRECTORY TO FILE LOCATION
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

#LOAD DATA
paper_risp_f<-read.dta("immigration_speeders.dta")

#EXPERIMENTAL FACTORS
paper_risp_f$approved_asy  <- factor(paper_risp_f$approved_asy)
paper_risp_f$exp_1  <- factor(paper_risp_f$exp_1)
paper_risp_f$exp_2  <- factor(paper_risp_f$exp_2)
paper_risp_f$exp_3  <- factor(paper_risp_f$exp_3)
paper_risp_f$ideo  <- factor(paper_risp_f$ideo)
paper_risp_f$ideo2  <- factor(paper_risp_f$ideo2)
paper_risp_f$Rotation_Block_F_IMM_T4  <- factor(paper_risp_f$Rotation_Block_F_IMM_T4)

#GLM 2  MAIN EFFECTS BY POSITION (NOT WEIGHTED)
m2way <- glm(approved_asy ~ exp_1*Rotation_Block_F_IMM_T4+exp_2*Rotation_Block_F_IMM_T4+exp_3*Rotation_Block_F_IMM_T4, family=binomial(), data = paper_risp_f)
summary(m2way)

#PREDICTIONS 2WAY
prediction(m2way, at = list(Rotation_Block_F_IMM_T4 = c("1st position", "2nd position", "3rd position"), exp_1 = c("Low skilled", "Skilled"))) %>% 
  summary -> pred2wayskills

prediction(m2way, at = list(Rotation_Block_F_IMM_T4 = c("1st position", "2nd position", "3rd position"), exp_2 = c("Looking for job", "Fleeing from war"))) %>% 
  summary -> pred2wayreason

prediction(m2way, at = list(Rotation_Block_F_IMM_T4 = c("1st position", "2nd position", "3rd position"), exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> pred2wayreasonorigin

### GGPLOT 2WAY
ggplot(data = pred2wayskills,
       aes(x = `at(exp_1)`,
           y = Prediction,
           shape = `at(Rotation_Block_F_IMM_T4)`,
           color = `at(Rotation_Block_F_IMM_T4)`)) +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0.5)
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue", "seagreen4")) +
  scale_shape_manual(name="", values = c(16, 17, 15)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(legend.title = element_text(size = 14),
        legend.text = element_text(size = 14)
  ) +
  theme(strip.text = element_text(size=14)) -> plot1

ggplot(data = pred2wayreason,
       aes(x = `at(exp_2)`,
           y = Prediction,
           shape = `at(Rotation_Block_F_IMM_T4)`,
           color = `at(Rotation_Block_F_IMM_T4)`)) +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0.5)
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue", "seagreen4")) +
  scale_shape_manual(name="", values = c(16, 17, 15)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(legend.title = element_text(size = 14),
        legend.text = element_text(size = 14)
  ) +
  theme(strip.text = element_text(size=14)) -> plot2

ggplot(data = pred2wayreasonorigin,
       aes(x = `at(exp_3)`,
           y = Prediction,
           shape = `at(Rotation_Block_F_IMM_T4)`,
           color = `at(Rotation_Block_F_IMM_T4)`)) +
  geom_point(size = 3, position = position_dodge(width = 0.5)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0.5)
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue", "seagreen4")) +
  scale_shape_manual(name="", values = c(16, 17, 15)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(legend.title = element_text(size = 14),
        legend.text = element_text(size = 14)
  ) +
  theme(strip.text = element_text(size=14)) -> plot3

g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

mylegend<-g_legend(plot1)
plotcombined <- grid.arrange(arrangeGrob(plot1 + theme(legend.position="none"),
                                         plot2 + theme(legend.position="none"),
                                         plot3 + theme(legend.position="none"),
                                         nrow=1),
                             mylegend, nrow=2,heights=c(10, 1))

ggsave("figureApp5.png", dpi = 600, plotcombined)


#GLM MAIN EFFECTS CONDITIONING BY THE ORDER (NOT WEIGHTED)
m3way <- glm(approved_asy ~ exp_1*exp_2*exp_3+Rotation_Block_F_IMM_T4, family=binomial(), data = paper_risp_f)
summary(m3way)

###PREDICTIONS MAIN EFFECTS
prediction(m3way, at = list(exp_1 = c("Low skilled", "Skilled"))) %>% 
  summary -> main_skilled

prediction(m3way, at = list(exp_2 = c("Looking for job", "Fleeing from war"))) %>% 
  summary -> main_reason

prediction(m3way, at = list(exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> main_origin


### GGPLOT MAIN EFFECTS CONDITIONING BY THE ORDER
ggplot(data = main_skilled,
       aes(x = `at(exp_1)`,
           y = Prediction,
           shape = `at(exp_1)`,
           color = `at(exp_1)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ylim(0, 1) +
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot1

ggplot(data = main_reason,
       aes(x = `at(exp_2)`,
           y = Prediction,
           shape = `at(exp_2)`,
           color = `at(exp_2)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("seagreen4", "palevioletred4")) +
  scale_shape_manual(name="", values = c(15, 18)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot2

ggplot(data = main_origin,
       aes(x = `at(exp_3)`,
           y = Prediction,
           shape = `at(exp_3)`,
           color = `at(exp_3)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("steelblue4", "antiquewhite4")) +
  scale_shape_manual(name="", values = c(4, 8)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot3

plot1 + plot2 + plot3

ggsave("figureApp6.png", dpi = 600)


#GLM MAIN EFFECTS (WEIGHTED)
m3wayw <- glm(approved_asy ~ exp_1*exp_2*exp_3, family = quasibinomial, weight = int_wt1_cap_T4, data = paper_risp_f)
summary(m3wayw)

###PREDICTIONS MAIN EFFECTS
prediction(m3way, at = list(exp_1 = c("Low skilled", "Skilled"))) %>% 
  summary -> main_skilled

prediction(m3way, at = list(exp_2 = c("Looking for job", "Fleeing from war"))) %>% 
  summary -> main_reason

prediction(m3way, at = list(exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> main_origin


### GGPLOT MAIN EFFECTS 
ggplot(data = main_skilled,
       aes(x = `at(exp_1)`,
           y = Prediction,
           shape = `at(exp_1)`,
           color = `at(exp_1)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot1

ggplot(data = main_reason,
       aes(x = `at(exp_2)`,
           y = Prediction,
           shape = `at(exp_2)`,
           color = `at(exp_2)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("seagreen4", "palevioletred4")) +
  scale_shape_manual(name="", values = c(15, 18)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot2

ggplot(data = main_origin,
       aes(x = `at(exp_3)`,
           y = Prediction,
           shape = `at(exp_3)`,
           color = `at(exp_3)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("steelblue4", "antiquewhite4")) +
  scale_shape_manual(name="", values = c(4, 8)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot3

plot1 + plot2 + plot3 

ggsave("figureApp7.png", dpi = 600)


#####MODELS WITHOUT SPEEDERS#####

#LOAD DATA
paper_risp_f<-read.dta("immigration_nospeeders.dta")

#EXPERIMENTAL FACTORS
paper_risp_f$approved_asy  <- factor(paper_risp_f$approved_asy)
paper_risp_f$exp_1  <- factor(paper_risp_f$exp_1)
paper_risp_f$exp_2  <- factor(paper_risp_f$exp_2)
paper_risp_f$exp_3  <- factor(paper_risp_f$exp_3)
paper_risp_f$ideo  <- factor(paper_risp_f$ideo)
paper_risp_f$ideo2  <- factor(paper_risp_f$ideo2)

#GLM MAIN EFFECTS (NOT WEIGHTED)
m3way <- glm(approved_asy ~ exp_1*exp_2*exp_3, family=binomial(), data = paper_risp_f)
summary(m3way)

###PREDICTIONS MAIN EFFECTS
prediction(m3way, at = list(exp_1 = c("Low skilled", "Skilled"))) %>% 
  summary -> main_skilled

prediction(m3way, at = list(exp_2 = c("Looking for job", "Fleeing from war"))) %>% 
  summary -> main_reason

prediction(m3way, at = list(exp_3 = c("Syrian", "Ukrainian"))) %>% 
  summary -> main_origin

### GGPLOT MAIN EFFECTS 
ggplot(data = main_skilled,
       aes(x = `at(exp_1)`,
           y = Prediction,
           shape = `at(exp_1)`,
           color = `at(exp_1)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("firebrick2", "cornflowerblue")) +
  scale_shape_manual(name="", values = c(16, 17)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot1

ggplot(data = main_reason,
       aes(x = `at(exp_2)`,
           y = Prediction,
           shape = `at(exp_2)`,
           color = `at(exp_2)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("seagreen4", "palevioletred4")) +
  scale_shape_manual(name="", values = c(15, 18)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot2

ggplot(data = main_origin,
       aes(x = `at(exp_3)`,
           y = Prediction,
           shape = `at(exp_3)`,
           color = `at(exp_3)`)) +
  geom_point(size = 3, position = position_dodge(width = 0)) +
  geom_pointrange(
    width = .1,
    aes(ymin = lower, ymax = upper),
    position = position_dodge(width = 0) 
  ) +
  #ggsci::scale_color_startrek(name = "") +
  theme_bw() +
  ylab("Asylum application should be approved (probability)") +
  xlab("") +
  theme(legend.position="bottom") +
  scale_colour_manual(name="", values = c("steelblue4", "antiquewhite4")) +
  scale_shape_manual(name="", values = c(4, 8)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + 
  ylim(0, 1) + 
  theme(axis.text=element_text(size=14),
        axis.title=element_text(size=14,face="bold")) +
  theme(strip.text = element_text(size=14)) + 
  theme(legend.position = "none") -> plot3

plot1 + plot2 + plot3 

ggsave("figureApp8.png", dpi = 600)
